! Fluxes Calculation
subroutine CalcFluxes()
use Data_Mesh, only: nCells, nFaces, meshFaceOwner, meshFaceNeighbour, nInternalFaces, &
                     nBoundaryFields, boundaryMeshName, boundaryMeshStart, boundaryMeshSize, &
                     df,db, meshCells
use DataCells, only: CharsCellPresent, SubStepCellVars
use DataFaces, only: CharsFacePresent, fluxesFacePresent, CharsNew, CharsOld, Acoef, &
                     fluxFacesReserve, fluxFaces, fluxCorrection, FaceVarsNew
use Data_Cfoam, only: btype
use DataPerfectGas, only: PressureAllVars, SoundVelocityPRho, invGammX, DensityPS, EnergyPRhoU
use DataTime, only: TimeCell, dtCell, TimeEpsilon, dtFace, dt2Cell, CurrentTimeCell, &
                    NextTimeCell, NextTimeLocal, nNextTimeCell,nNextTimeLocal, &
                    nNextTimeCell1
use DataViscosity, only: CalcViscosityOn
implicit none

    integer :: i, iCell, iBoundField, iFace, iCellOwner,iCellNeighbour
    logical :: CalcAcoef, CalcNewCharsOn(2)
    real(8) :: NoSlipCoef
    real(8) :: CharsLocal(5,2), CharsLocalNew(5)
    real(8) :: EigenVal3SCF,EigenVal3SCB, SoundSCF,SoundSCB
    real(8) :: AcoefQNew,AcoefRNew
    real(8) :: Time10,Time20, Time1,Time2, dt(2), Time, dt1,dt2
    real(8) :: ds1,ds2
    real(8) :: Fluxes(5)
    real(8) :: U(3), Pres, Rho

!        if(NextTimeFace(iFace).gt.NextTimeLocal+TimeEpsilon) cycle

!    do iFace = 1,nInternalFaces

!    NextTimeLocal = minVal(NextTimeCell(1:nCells))
    nNextTimeLocal = minVal(nNextTimeCell(1:nCells))
    !SK nNextTimeCell1(1:nCells) = nNextTimeCell(1:nCells) - nNextTimeLocal ! ==0 = .false./ !=0 = .true.
    !do iCell = 1,nCells
     !   nNextTimeCell1(iCell) = nNextTimeCell(iCell) - nNextTimeLocal ! ==0 = .false./ !=0 = .true.
    !end do

     do iCell = 1,nCells
        if(nNextTimeCell(iCell) - nNextTimeLocal.eq.0)then
            nNextTimeCell1(iCell) = .false. ! ==0 = .false./ !=0 = .true.
        else
            nNextTimeCell1(iCell) = .true. ! ==0 = .false./ !=0 = .true.
        end if
    end do


!    if(nNextTimeLocal.lt.0)then
!        write(*,*) '@@@ nNextTimeLocal=', nNextTimeLocal
!        Stop 'nNextTimeLocal<0'
!    end if

    do iCell = 1,nCells
!        if(NextTimeCell(iCell).gt.NextTimeLocal+TimeEpsilon) cycle
!        if(nNextTimeCell(iCell).ne.nNextTimeLocal) cycle

        if(nNextTimeCell1(iCell)) cycle
        
        do i = 1,6
            iFace = meshCells(i,iCell)
            if(iFace.le.nInternalFaces)then
                if(fluxesFacePresent(1,iFace).or.fluxesFacePresent(2,iFace)) cycle
            
                iCellNeighbour = meshFaceNeighbour(iFace)
                iCellOwner = meshFaceOwner(iFace)

                EigenVal3SCB = 0.5d0 * (CharsNew(1,1,iFace)+CharsNew(2,1,iFace))
                EigenVal3SCF = 0.5d0 * (CharsNew(1,2,iFace)+CharsNew(2,2,iFace))

                Pres = PressureAllVars(SubStepCellVars(1:5,iCellOwner))
                Rho = SubStepCellVars(1,iCellOwner)
                SoundSCB = SoundVelocityPRho(Pres,Rho)

                Pres = PressureAllVars(SubStepCellVars(1:5,iCellNeighbour))
                Rho = SubStepCellVars(1,iCellNeighbour)
                SoundSCF = SoundVelocityPRho(Pres,Rho)

                CharsLocal(1:5,1:2) = CharsNew(1:5,1:2,iFace)

                ds1 = db(iFace)
                ds2 = df(iFace)

                ! Definition Interpolated Characteristis on the Face
                Time10 = TimeCell(iCellOwner)
                Time20 = TimeCell(iCellNeighbour)

                Time1 = TimeCell(iCellOwner) + dtCell(iCellOwner)
                Time2 = TimeCell(iCellNeighbour) + dtCell(iCellNeighbour)

                dt(1) = dtCell(iCellOwner)
                dt(2) = dtCell(iCellNeighbour)
                CalcNewCharsOn(1:2) = .false.
        !        fluxesFacePresent(1:2,iFace) = .true.

                if(abs(Time1-Time2).lt.TimeEpsilon)then
                    fluxesFacePresent(1:2,iFace) = .true.
                    CalcNewCharsOn(1:2) = .true.
                    CharsLocal(1:5,1:2) = CharsNew(1:5,1:2,iFace)
                else
                    if(Time1.lt.Time2)then
                        CharsLocal(1:5,1) = CharsNew(1:5,1,iFace)
                        call CharsLinearInterpolation(Time20,Time1,Time2,2)
                        fluxesFacePresent(1,iFace) = .true.
                        CalcNewCharsOn(1) = .true.
                    end if

                    if(Time2.lt.Time1)then
                        CharsLocal(1:5,2) = CharsNew(1:5,2,iFace)
                        call CharsLinearInterpolation(Time10,Time2,Time1,1)
                        fluxesFacePresent(2,iFace) = .true.
                        CalcNewCharsOn(2) = .true.
                    end if
                end if

                call CalcFaceFluxes(EigenVal3SCF,SoundSCF, EigenVal3SCB,SoundSCB)

                ! Definition of Conservative Flux Correction
        !        if(CalcNewCharsOn(1))then
                if(fluxesFacePresent(1,iFace)) then
                    dt1 = max(0d0,min(dt2Cell(iCellOwner),Time10 + dt2Cell(iCellOwner) - Time20))
                    dt2 = min(dt2Cell(iCellOwner),Time1 - Time20)
                    fluxFacesReserve(1:5,2,iFace) = fluxFacesReserve(1:5,2,iFace) + &
                                               dt1 * (-fluxFaces(1:5,1,iFace)) + &
                                               dt2 * (-Fluxes(1:5))
                end if

        !        if(CalcNewCharsOn(2))then
                if(fluxesFacePresent(2,iFace)) then
                    dt1 = max(0d0,min(dt2Cell(iCellNeighbour),Time20 + dt2Cell(iCellNeighbour) - Time10))
                    dt2 = min(dt2Cell(iCellNeighbour),Time2 - Time10)
                    fluxFacesReserve(1:5,1,iFace) = fluxFacesReserve(1:5,1,iFace) + &
                                               dt1 * (-fluxFaces(1:5,2,iFace)) + &
                                               dt2 * (Fluxes(1:5))
                end if

        !        if(CalcNewCharsOn(1))then
                if(fluxesFacePresent(1,iFace)) then
        !            if(TimeCell(iCellOwner).lt.TimeCell(iCellNeighbour))then
                    if(TimeCell(iCellNeighbour)-TimeCell(iCellOwner).gt.TimeEpsilon)then
                        if(fluxesFacePresent(2,iFace)) then
                            dt1 = dt2Cell(iCellOwner)
                            dt2 = dt1

                            dt1 = min(dt2Cell(iCellOwner),max(0d0,Time2-Time10))
                            dt2 = max(0d0,Time2-Time10-dt2Cell(iCellOwner))
                        else
                            dt1 = min(dt2Cell(iCellOwner),max(0d0,Time20-Time10))
                            dt2 = max(0d0,Time20-Time10-dt2Cell(iCellOwner))
                        end if
                        if(dt1+dt2.gt.TimeEpsilon) then
                            fluxCorrection(1:5,1,iFace) = (fluxFacesReserve(1:5,1,iFace) - &
                                                 (dt1 * fluxFaces(1:5,1,iFace) + dt2 * Fluxes(1:5)) )
                        end if
                    end if
                    fluxFacesReserve(1:5,1,iFace) = 0d0

                    fluxFaces(1:5,1,iFace) = Fluxes(1:5)
                    FaceVarsNew(1,1,iFace) = Rho
                    FaceVarsNew(2:4,1,iFace) = U(1:3)
                    FaceVarsNew(5,1,iFace) = Pres ! Energy
                end if

        !        if(CalcNewCharsOn(2))then
                if(fluxesFacePresent(2,iFace)) then
        !            if(TimeCell(iCellOwner).gt.TimeCell(iCellNeighbour))then
                    if(TimeCell(iCellOwner)-TimeCell(iCellNeighbour).gt.TimeEpsilon)then
                        if(fluxesFacePresent(1,iFace)) then
                            dt1 = dt2Cell(iCellNeighbour)
                            dt2 = dt1

                            dt1 = min(dt2Cell(iCellNeighbour),max(0d0,Time1-Time20))
                            dt2 = max(0d0,Time1-Time20-dt2Cell(iCellNeighbour))
                        else
                            dt1 = min(dt2Cell(iCellNeighbour),max(0d0,Time10-Time20))
                            dt2 = max(0d0,Time10-Time20-dt2Cell(iCellNeighbour))
                        end if
                        if(dt1+dt2.gt.TimeEpsilon) then
                            fluxCorrection(1:5,2,iFace) = (fluxFacesReserve(1:5,2,iFace) - &
                                                 (dt1 * fluxFaces(1:5,2,iFace) - dt2 * Fluxes(1:5)) )
                        end if
                    end if
                    fluxFacesReserve(1:5,2,iFace) = 0d0
                    fluxFaces(1:5,2,iFace) = -Fluxes(1:5)
                    FaceVarsNew(1,2,iFace) = Rho
                    FaceVarsNew(2:4,2,iFace) = U(1:3)
                    FaceVarsNew(5,2,iFace) = Pres ! Energy
                end if
!            end do
            
            else

!    do iFace = nInternalFaces+1,nFaces
                if(fluxesFacePresent(1,iFace)) cycle

                iCellOwner = meshFaceOwner(iFace)

                EigenVal3SCF = 0.5d0 * (CharsNew(1,2,iFace)+CharsNew(2,2,iFace))
                EigenVal3SCB = 0.5d0 * (CharsNew(1,1,iFace)+CharsNew(2,1,iFace))

                Pres = PressureAllVars(SubStepCellVars(1:5,iCellOwner))
                Rho = SubStepCellVars(1,iCellOwner)
                SoundSCB = SoundVelocityPRho(Pres,Rho)

                Pres = FaceVarsNew(5,2,iFace)
                Rho = FaceVarsNew(1,2,iFace)
                SoundSCF = SoundVelocityPRho(Pres,Rho)

                CharsLocal(1:5,1:2) = CharsNew(1:5,1:2,iFace)

                ds1 = db(iFace)
                ds2 = ds1
                
                CalcNewCharsOn(1:2) = .true.

                call CalcFaceFluxes(EigenVal3SCF,SoundSCF, EigenVal3SCB,SoundSCB)

                fluxFaces(1:5,1,iFace) = Fluxes(1:5)
                fluxesFacePresent(1,iFace) = .true.
                FaceVarsNew(1,1,iFace) = Rho
                FaceVarsNew(2:4,1,iFace) = U(1:3)
                FaceVarsNew(5,1,iFace) = Pres ! Energy
                
                FaceVarsNew(1:5,2,iFace) = FaceVarsNew(1:5,1,iFace)
!        end do
            end if
        end do
    end do
    
contains

    ! Selection which R and Q characteristics should be taken (left or right)
    subroutine CalcFaceFluxes(EigenVal3SCF,SoundSCF, EigenVal3SCB,SoundSCB)
    implicit none
    real(8),intent(in) :: EigenVal3SCF,SoundSCF, EigenVal3SCB,SoundSCB

        real(8) :: EigenVal1SCF,EigenVal1SCB, EigenVal2SCF,EigenVal2SCB

        EigenVal1SCF = EigenVal3SCF - SoundSCF
        EigenVal1SCB = EigenVal3SCB - SoundSCB
        call CalcCharsNewMinVar(EigenVal1SCF,EigenVal1SCB,AcoefQNew,1)
        ! SK END OF MINIMUM VAR CHAR1

        ! SK MINIMUM VARIATION -- CHAR2
        EigenVal2SCF = EigenVal3SCF + SoundSCF
        EigenVal2SCB = EigenVal3SCB + SoundSCB
        call CalcCharsNewMinVar(EigenVal2SCF,EigenVal2SCB,AcoefRNew,2)
        ! SK END OF MINIMUM VAR CHAR2

        call CalcLocalFluxes()

    end subroutine CalcFaceFluxes

    ! Selection which one characteristic should be taken
    subroutine CalcCharsNewMinVar(EigenValSCF,EigenValSCB,AcoefNew,iChar)
    use Data_cfoam, only: Usnew, Psnew, Rhosnew, Tsnew
    use Data_Mesh, only: meshFaceCenters
    use DataTime, only: nStepLocal,TimeLocal, nNextTimeLocal, ttime
    implicit none
        real(8),intent(in) :: EigenValSCF,EigenValSCB
        integer,intent(in) :: iChar
        real(8),intent(out) :: AcoefNew

        if(EigenValSCF.gt.0d0.and.EigenValSCB.gt.0d0)then
            CharsLocalNew(iChar) = CharsLocal(iChar,1)
            AcoefNew = Acoef(1,iFace) ! AcoefNew = AcoefB
            return
        end if

        if(EigenValSCF.lt.0d0.and.EigenValSCB.lt.0d0)then
            CharsLocalNew(iChar) = CharsLocal(iChar,2)
            AcoefNew = Acoef(2,iFace) ! AcoefNew = AcoefF
            return
        end if

        write(*,*) 'meshFaceCenters', meshFaceCenters(1:3,iFace)
        write(*,*) 'Vars1=', FaceVarsNew(1:5,1,iFace)
        write(*,*) 'Vars2=', FaceVarsNew(1:5,2,iFace)
        write(*,'(i8,2e20.10,2i12)') nStepLocal,ttime+TimeLocal,minVal(dtCell(:)), nNextTimeLocal, iFace
        STOP '@@ Note: Slipping place! @@'

    end subroutine CalcCharsNewMinVar

    ! Calcualtion of fluxes through face iFace
    subroutine CalcLocalFluxes()
    use Data_Mesh, only: meshSf
    implicit none

        real(8) :: Usnew1D, Ut(3),UtNew(3), alfa, gam, a(3,3), n(3),normal(3)
        integer :: iSide
        real(8) :: Un, Energy, Umod, sound 

        Usnew1D = CharsLocalNew(1) + AcoefQNew*(CharsLocalNew(2)-CharsLocalNew(1))/(AcoefRNew + AcoefQNew)
        !qn+gf*(rn-qn)/(gf+gb)

        call GetAlfa(Usnew1D,alfa)

        if(alfa>=1.0d0) then   ! if(Usnew1D<0.0) then ! SK changed, 2011-05-17
            CharsLocalNew(3:5) = CharsLocal(3:5,2)
        else if(alfa<=0.0d0) then
            CharsLocalNew(3:5) = CharsLocal(3:5,1) ! if(Usnew1D>0.0)
        else
            CharsLocalNew(3:5) = alfa * CharsLocal(3:5,2) + (1d0-alfa) * CharsLocal(3:5,1)
        end if

        do iSide = 1,2
            if(CalcNewCharsOn(iSide).eqv..false.) cycle
!            CharsOld(1:5,iSide,iFace) = CharsNew(1:5,iSide,iFace)
            CharsNew(1:5,iSide,iFace) = CharsLocalNew(1:5)
        end do

        ! Definition of Variables on Face on New Time Level
        Pres = ((CharsLocalNew(2)-CharsLocalNew(1)) / (AcoefRNew + AcoefQNew))**invGammX
!        Pres = (0.5d0 * (CharsLocalNew(2)-CharsLocalNew(1)) / AcoefRNew )**invGammX
        Rho = DensityPS(Pres,CharsLocalNew(3))

        Ut(1) = 0d0
        Ut(2:3) = CharsLocalNew(4:5)

        n(1:3) = meshSf(1:3,iFace)
        normal(1:3) = n(1:3)/dsqrt(sum(n(1:3)*n(1:3)))
        ! transformation {x,y,z} -> {normal,tan1,tan2}
        call coor(normal(1),normal(2),normal(3),a)
        call transformr(Ut,UtNew,a)

        U(1:3)=UtNew(1:3)+Usnew1D*normal(1:3)
        
        Energy = EnergyPRhoU(Pres,Rho,U)

	    Umod = dsqrt(sum(U(1:3)*U(1:3)))
        sound = SoundVelocityPRho(Pres,Rho)
        dtFace(iFace) = min(abs(ds1),abs(ds2)) / (Umod + sound) ! NOTE: This is needed to be take into account!
!	    if (dtmint.lt.dtmin) dtmin=dtmint

        ! Definition of Fluxes on New Time Level
        Un = sum(n(1:3)*U(1:3))
        Fluxes(1) = Un * Rho
        Fluxes(2:4) = Un*Rho*U(1:3) + (Pres-1d5)*n(1:3) ! Pres*n(1:3) ! 
!        Fluxes(2:4) = Un*Rho*U(1:3) + Pres*n(1:3) ! Pres*n(1:3) ! 
        Fluxes(5) = Un * (Rho*Energy + Pres)

    end subroutine CalcLocalFluxes

    ! Interpolation of charachteristics
    subroutine CharsLinearInterpolation(Time0,TimeMiddle,Time1,iSide)
    implicit none
        real(8),intent(in) :: Time0,TimeMiddle,Time1
        integer,intent(in) :: iSide

        CharsLocal(1:5,iSide) = CharsOld(1:5,iSide,iFace) + (TimeMiddle-Time0) / (Time1 - Time0) * &
                                (CharsNew(1:5,iSide,iFace) - CharsOld(1:5,iSide,iFace))

    end subroutine CharsLinearInterpolation

end subroutine CalcFluxes

